Date read in

We used log density to standardize the variance in protists density time series. Since in some days, the protists density is 0, and will result -INF by log(), we added 1 to all the density before calculating the log density for model fitting. We minused the 1 after exponentiating the model predictions.

### data read in------
## density 
milou_dens0 <- read.csv("/Users/zeyihan/Documents/PhD/Projects/Milou/02_Clean_data/Milou_density.csv")
milou_dens <- milou_dens0  %>% mutate(tre=as.factor(tre), temp = as.factor(temp), nut = as.factor(nut), all_rep = rep(1:24, time = 24), log_density = log(density + 1)) %>% select(-1)  
tetra <- milou_dens[which(milou_dens$Class=="Tetrahymena"),]
eup <- milou_dens[which(milou_dens$Class=="Euplotes"),]
bac0 <- read.csv("/Users/zeyihan/Documents/PhD/Projects/Milou/02_Clean_data/Milou_bac_OD.csv")
bac1 <- bac0 %>% mutate(tre = as.factor(tre), temp = as.factor(temp), nut = as.factor(nut), all_rep = rep(1:24, time = 11), temp_num = as.numeric(rep(c(22,25), each = 6, times = 22)), nut_num = as.numeric(rep(c(1,0.5), each = 12, times = 11)))

We calculated the mean and standard variation for bacteria OD after day 5. One data point that is more than 3 standard deviation away from the mean is removed.

sd <- sd(bac1$OD[which(bac1$day > 5)], na.rm = TRUE)
mean <- mean(bac1$OD[which(bac1$day > 5)], na.rm = TRUE)
bac <- bac1[-which(bac1$OD > mean+3*sd),]

Because after population collapse, Tetrahymena are functionally extinct and so, even if we occasionally found individual, trait distributions are unreliable. For traits observation made after each Tetrahymena treatment collapsed (once population density became 0 ind/ml) that with less than 10 number of individual, we considered them unreliable and, thus, removed for future analysis. For time series analysis, we only used Tetrahymena trait data from day 0-9 since the analysis we used required equal-length time series.

Figure 1

Calculate mean density of Bacteria, Tetrahymena, and Euplotes

bac.2 <- bac %>%
  rename(density = OD) %>% 
  select(tre:density,treatNum) %>%
  mutate(temp = as.factor(temp))

teb.all <- milou_dens %>% 
 select(tre:density,temp_num,-temp,treatNum)  %>%
  rename(temp = temp_num) %>%
  mutate(temp = as.factor(temp), nut = as.factor(nut),day = as.numeric(day))%>%
  full_join(bac.2) %>% 
  na.omit() %>%
  group_by(Class,tre,treatNum, nut, temp,day) %>%
  summarise(mean.d = mean(density),
            log.d = log(mean.d +1)) %>%
  mutate(day = as.numeric(day))

teb.mean <- teb.all %>%
 select(Class:mean.d) %>%
  spread(key=Class, value = mean.d)

Fig 1a. Overall average ecological dynamics

dens1 <- tetra %>% group_by(day) %>% summarise(tetra = log(mean(density+1)), t.err = sd(density+1)/sqrt(24), t.err.log = log(t.err))
d.eup <- eup %>% group_by(day) %>% summarise(eup = log(mean(density+1)), eup.err = sd(density+1)/sqrt(24), eup.err.log = log(eup.err) ) 
dens2 <- right_join(dens1,d.eup)
d.bac <- bac %>% group_by(day) %>% na.omit() %>% summarise(bac = mean(OD), bac.err = sd(OD)/sqrt(24)) 
dens <- full_join(dens2,d.bac)
# Tetrahymena and Euplotes average abundance
par(mar = c(5,5,2,5),tcl = 0.5)
with(dens, plot(day, tetra, type="o",pch=3, cex =2, lwd = 4,col="#818181", ylab="Tetrahymena/Euplotes Abundance", ylim=c(0,12), axes = F))
axis(side = 1,las=1,cex.axis=1.5,lwd = 2)
axis(side = 2,las=1,cex.axis=1.5,lwd = 2)
lines(dens$day,dens$eup, col = "black", lwd = 4) #lty = 2
points(dens$day,dens$eup, col = "black", pch =4, cex =2,lwd = 4)
# Bacteria average OD
par(new = T,tcl = -0.5)
with(dens, plot(day,bac, type="o", col ="#bcbcbc", cex=2, ylim=c(0, 0.03), lwd = 4, pch = 8,axes=F, xlab=NA, ylab=NA)) 
axis(side = 4, at = c(0, 0.01,0.02, 0.03),las=1,cex.axis=1.5,lwd = 2)
mtext(side = 4, 'Bacteria OD')

Fig 1b-d. Gamm fits for the ecological dynamics

Based on model comparison (Supplementray Table X), we used the Generalized Additive Mixed Models (GAMM) with 4 treatments for all species for both ecological and phenotypic dynamics.

# Bacteria ecological dyanmics
bm2 <- gamm(OD ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = bac, correlation = corARMA(form = ~ day|tre/rep, p = 2)) 
# Tetrahymena ecological dynamics
tm1 <- gamm(log_density ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = tetra, correlation = corARMA(form = ~ day|tre/rep, p = 1)) 
# Euplotes ecological dynamics
em0 <- gamm(log_density ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = eup) 

# function to predict tetrahymena, euplotes, and bacteria density from previous selected models
em0 <- gamm(log_density ~ s(day, by = tre, k =10) + tre,random = list(rep= ~ 1), data = eup) 
tm1 <- gamm(log_density ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = tetra, correlation = corARMA(form = ~ day|tre/rep, p = 1)) 
bm2 <- gamm(OD ~ s(day, by = tre, k =10) + tre ,random = list(rep= ~ 1), data = bac, correlation = corARMA(form = ~ day|tre/rep, p = 2)) 

pred_gam_TEB <- function(treatment,day, length){
  df <- data.frame(day=seq(0, day,length.out=length), tre=rep(treatment,length))
  eup_pred <-  exp(predict.gam(em0$gam,df,type = "response", se.fit =T)$fit)-1
  tetra_pred <-  exp(predict.gam(tm1$gam,df,type = "response", se.fit =T)$fit)-1 
  bac_pred <- predict.gam(bm2$gam,df,type = "response", se.fit =T)$fit
  pred <- as.data.frame(cbind(df[, 1], eup_pred, tetra_pred, bac_pred))
  colnames(pred) <- c("day","eup_pred", "tetra_pred", "bac_pred" )
  pred
}

F22gam <- pred_gam_TEB("Full22C",15, 800)
F25gam <- pred_gam_TEB("Full25C",15, 800)
H22gam <- pred_gam_TEB("Half22C",15, 800)
H25gam <- pred_gam_TEB("Half25C",15, 800)

gam_all <- rbind(F22gam,F25gam,H22gam,H25gam) %>% 
  mutate(tre=rep(c("Full22C", "Full25C", "Half22C", "Half25C"),each = 800),
         temp = rep(c(22,25,22,25), each=800),
         nut = rep(c("Full", "Half"), each=1600),
         day=rep(seq(0,15,length.out=800),4)) 
gam_gather <- gam_all%>%
  gather(key = "Class", value = "density", eup_pred, tetra_pred , bac_pred)

Fig 1b Bacteria overall ecological dynamics

trt_name <- c(expression(paste("Full 22 ",degree,"C")),expression(paste("Half 22 ",degree,"C")), expression(paste("Full 25 ",degree,"C")),expression(paste("Half 25 ",degree,"C")) )

# Bacteria
ggplot(data = bac, aes(x=day,y=OD,group = temp,shape = tre,colour = tre))+ 
  geom_line(data = gam_all, aes(x=day,y=bac_pred,group = tre, colour = tre), size = 0.8)+ 
  geom_jitter(size = 1.5, width = 0.7)+
  scale_colour_manual(name = "Treatment",labels = trt_name,values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment", labels = trt_name, values = c(19, 19, 17, 17))+
  theme_bw()+theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.ticks.length=unit(-0.2, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")))

Fig 1c Tetrahymena overall ecological dynamics

# Tetrahymena
ggplot(data = tetra, aes(x=day,y=density,group = temp,shape = tre,colour = tre))+ 
  geom_line(data = gam_all, aes(x=day,y=tetra_pred,group = tre, colour = tre), size = 0.8)+  
  geom_jitter(size = 1.5, width = 0.7)+
  scale_colour_manual(name = "Treatment",labels = trt_name,values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment", labels = trt_name, values = c(19, 19, 17, 17))+
  theme_bw()+theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.ticks.length=unit(-0.2, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")))

Fig 1d Euplotes overall ecological dynamics

# Euplotes
ggplot(data = eup, aes(x=day,y=density,group = temp,shape = tre,colour = tre)) +
  geom_line(data = gam_all, aes(x=day,y=eup_pred,group = tre, colour = tre), size = 0.8)+ 
geom_jitter(size = 1.5, width = 0.7)+
  scale_colour_manual(name = "Treatment",labels = trt_name,values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment", labels = trt_name, values = c(19, 19, 17, 17))+
  theme_bw()+theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.ticks.length=unit(-0.2, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")))

Fig 1e-f GAMM fits for phenotypic dynamics of Tetrahymena and Euplotes

Fig 1e Tetrahymena overall phenotypic dynamics

tt3 <- gamm(area ~ s(day, by = tre, k =7) + tre,random = list(rep= ~ 1), data = tetraSumm, correlation = corARMA(form = ~ day|tre/rep, p = 3))
ttgam <- data.frame(day=rep(seq(0,9,length.out=500), 4), tre=rep(c("Full22C", "Full25C", "Half22C", "Half25C"),each = 500, 4))  
ttgam$tetra <-  predict.gam(tt3$gam,ttgam,type = "response", se.fit =T)$fit

ggplot(ttgam, aes(x=day, y=tetra, color = tre, shape = tre))+ geom_line(width = 2) +
  geom_jitter(data= tetraSumm, aes(x=day, y = area, color=tre),size = 1.3, jitter = 0.7)+
  #geom_vline(xintercept=9,color = "#000000",  size = 0.3)+
   scale_colour_manual(name = "Treatment",labels = trt_name,values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment", labels = trt_name, values = c(19, 19, 17, 17))+
  theme_bw()+theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.ticks.length=unit(-0.15, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")))

Fig 1f Euplotes overall phenotypic dynamics

# make days with number of individal less than 10 as NA
#eupSumm$area[which(eupSumm$num_ind < 11)] <- NA

et1 <- gamm(area ~ s(day, by=tre, k =9) + tre,random = list(rep= ~ 1), data = eupSumm, correlation = corARMA(form = ~ day|tre/rep, p = 1)) 
etgam <- data.frame(day=rep(seq(0,15,length.out=1000), 4), tre=rep(c("Full22C", "Full25C", "Half22C", "Half25C"),each = 1000, 4))  
etgam$eup<-  predict.gam(et1$gam,etgam,type = "response", se.fit =T)$fit

ggplot(etgam, aes(x=day, y=eup, color = tre, shape = tre))+ geom_line(width = 2) +
  geom_jitter(data= eupSumm, aes(x=day, y = area, color=tre),size = 1.3, jitter = 0.7)+
   scale_colour_manual(name = "Treatment",labels = trt_name,values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment", labels = trt_name, values = c(19, 19, 17, 17))+
  theme_bw()+theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.ticks.length=unit(-0.15, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")))

Fig 2 Interactive plot of temperature and nutrients effects on ecological dynamics

Descriptors for Species Ecological Dynamics

The following are ecological dynamics descriptor calculated for each species.

1. Bacteria ecological dynamics descriptor

  • Initial growth
  • Maximum abundance
  • Coefficient of variation (CV)
## New data frame
bac_desc <- data.frame("tre" = rep(c("Full22C", "Full25C", "Half22C", "Half25C"), each = 6))
bac_desc<- bac_desc %>%
  mutate(temp = as.factor(rep(c(22,25), each=6, times=2)), 
         nut = factor(rep(c("Full","Half"), each=12), levels = c("Half", "Full")), 
         rep = rep(1:6, times = 4))

## initial growth rate of Bacteria from day0 to day 1
bac_day0 <- bac[which(bac$day == 0), ]
bac_day1 <- bac[which(bac$day == 1), ]
bac_desc$init_growth <- (bac_day1$OD-bac_day0$OD)/8

## max abundance: three days with the top three population density
bac_max1<- teb.mean%>% group_by(tre,treatNum) %>% select(tre:day,Bacteria) %>% na.omit() %>%  arrange(desc(Bacteria), .by_group = TRUE)%>% slice_max(Bacteria, n=3) 

bm1 <- bac %>% filter(tre == "Full22C", day %in% bac_max1$day[which(bac_max1$tre=="Full22C")])
bm2 <- bac %>% filter(tre == "Full25C", day %in% bac_max1$day[which(bac_max1$tre=="Full25C")])
bm3 <- bac %>% filter(tre == "Half22C", day %in% bac_max1$day[which(bac_max1$tre=="Half22C")])
bm4 <- bac %>% filter(tre == "Half25C", day %in% bac_max1$day[which(bac_max1$tre=="Half25C")])

bac_max <- rbind(bm1, bm2, bm3,bm4) %>% select(tre:temp, OD, -Date) 
bac_max$nut <- factor(bac_max$nut, levels =  c("Half", "Full"))

# CV
bac_desc<- bac%>% group_by(tre,rep) %>%na.omit()%>% summarize(sd = sd(OD), meanOD= mean(OD), CV = sd/mean) %>% full_join(bac_desc) 

# mean
bac.mean_desc <- bac_desc %>% 
  group_by(tre, nut, temp) %>%
  na.omit %>%
  summarise(ig = mean(init_growth),
            CV.m=mean(CV))
bac_max_mean <- bac_max %>% group_by(nut, temp,tre) %>% na.omit() %>% summarise(ODmax = mean(OD),ODse = sd(OD))

The data frames for bacteria initial growth and coefficient of variation (CV) looks like:

head(bac_desc)
## # A tibble: 6 × 8
## # Groups:   tre [2]
##   tre       rep     sd meanOD    CV temp  nut   init_growth
##   <chr>   <int>  <dbl>  <dbl> <dbl> <fct> <fct>       <dbl>
## 1 Full22C     1 0.0105 0.0237 0.449 22    Full      0.00175
## 2 Full22C     2 0.0117 0.0251 0.500 22    Full      0.00213
## 3 Full22C     3 0.0126 0.0260 0.535 22    Full      0.00338
## 4 Full22C     4 0.0123 0.0270 0.525 22    Full      0.00275
## 5 Full22C     5 0.0109 0.0239 0.463 22    Full      0.00362
## 6 Full25C     1 0.0138 0.0302 0.586 25    Full      0.00312

We used the top three observed OD600 values for each replicate from each treatment for bacteria as their maximun abundance because their population dynamics remained relatively stable.The dataframe looks like:

head(bac_max_mean)
## # A tibble: 4 × 5
## # Groups:   nut, temp [4]
##   nut   temp  tre      ODmax    ODse
##   <fct> <fct> <fct>    <dbl>   <dbl>
## 1 Half  22    Half22C 0.0233 0.00475
## 2 Half  25    Half25C 0.0167 0.00153
## 3 Full  22    Full22C 0.0343 0.00689
## 4 Full  25    Full25C 0.0311 0.00951

2. Tetrahymena ecological dynamics descriptor

  • Initial growth
  • Maximum abundance
  • Coefficient of variation (CV)
  • Day of population collapse
# new data frame
tetra_desc <- data.frame("tre" = rep(c("Full22C", "Full25C", "Half22C", "Half25C"), each = 6))
tetra_desc<-tetra_desc %>%
  mutate(temp = as.factor(rep(c(22,25), each=6, times=2)), 
         nut = factor(rep(c("Full","Half"), each=12), levels = c("Half", "Full")), 
         rep = rep(1:6, times = 4))

## initial growth rate of tetrahymena from day0 to day 1 
tetra_day0 <- tetra%>% filter(day==0) %>% select(density)
tetra_day1 <- tetra%>% filter(day==1) %>% select(density)
tetra_desc$tetra.ig <- log(tetra_day1$density)-log(tetra_day0$density)

 ## max abundance
tetra_max1 <- teb.mean %>% group_by(tre,treatNum) %>% select(tre:day,Tetrahymena) %>% summarize(max = max(Tetrahymena))
tetra_max <- teb.mean %>% filter(Tetrahymena %in% tetra_max1$max)
tetra_desc <- tetra%>% filter(tre %in% tetra_max$tre & day %in% tetra_max$day)%>% select(tre,rep,density) %>% rename(max = density)%>% left_join(tetra_desc)

# CV
tetra_desc <- tetra%>% group_by(tre,rep) %>% summarize(sd = sd(density), mean= mean(density), CV = sd/mean)%>% full_join(tetra_desc)

## time collapse
tetra_desc <- tetra %>% group_by(tre,rep) %>% filter(density == 0) %>% summarize(tc = min(day)) %>% full_join(tetra_desc)

# mean of the descriptors
t.mean_desc <- tetra_desc %>%
  group_by(tre, nut, temp) %>%
  summarise(ig = mean(tetra.ig),
            max.m = mean(max), 
            CV.m = mean(CV),
            tc.m = mean(tc)) 

Data frame looks like:

head(tetra_desc)
## # A tibble: 6 × 10
## # Groups:   tre [1]
##   tre       rep    tc    sd  mean    CV   max temp  nut   tetra.ig
##   <chr>   <int> <int> <dbl> <dbl> <dbl> <dbl> <fct> <fct>    <dbl>
## 1 Full22C     1    10 2455. 1926.  1.27  6195 22    Full      2.77
## 2 Full22C     2    11 2841. 2309.  1.23  6354 22    Full      2.61
## 3 Full22C     3    14 2852. 2014.  1.42  7465 22    Full      3.12
## 4 Full22C     4    10 2599. 1911.  1.36  6966 22    Full      2.67
## 5 Full22C     5    10 2718. 2011.  1.35  6663 22    Full      2.82
## 6 Full22C     6    14 2704. 2545.  1.06  6733 22    Full      2.78

3. Euplotes ecological dynamics descriptor

  • Initial growth
  • Maximum abundance
  • Coefficient of variation (CV)
  • Day of population peak

Data frame looks like:

head(eup_desc)
## # A tibble: 6 × 11
## # Groups:   tre [1]
##   tre       rep    sd  mean    CV temp  nut   init_growth maxdp daypeak   max
##   <chr>   <int> <dbl> <dbl> <dbl> <fct> <fct>       <dbl> <dbl>   <int> <dbl>
## 1 Full22C     1  93.9  88.7  1.06 22    Full        0.782   246      10   227
## 2 Full22C     2  93.3  82.0  1.14 22    Full        0.675   250      14   250
## 3 Full22C     3 101.   87.6  1.15 22    Full        0.726   276      10   257
## 4 Full22C     4 104.   94.9  1.09 22    Full        0.738   284      14   284
## 5 Full22C     5  91.5  89.4  1.02 22    Full        0.815   270      14   270
## 6 Full22C     6 107.   89.6  1.19 22    Full        0.675   303      10   233

Fig 2 a-c Initial growth

Here is the function to generate interactive plot.

TEB_interactive <- function(data1,x1,y1,data2,x2,y2,y_label){
  p <- ggplot(data=data1, aes(x=x1, y=y1, group=temp, shape= tre,colour=tre)) +
    geom_line(data = data2, aes(x=x2, y=y2), size = 0.6, color = "#c3c3c3") +
    scale_x_discrete(limits=c("Half","Full")) +
    labs(y=y_label, x = "Nutrient") +
    geom_jitter(size = 2, shape = rep(c(1,2),each=12), width = 0.25) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width=0.2, size = 0.5) +
    geom_point(data = data2, aes(x=x2, y=y2, group=temp, shape=tre, colour=tre), size = 3) +
    theme_bw() +
    theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),axis.ticks.length=unit(-0.12, "cm"), axis.text.x = element_text(margin=unit(c(0.25,0.25,0.25,0.25), "cm"),colour = "black"), axis.text.y = element_text(margin=unit(c(0.25,0.25,0.25,0.25), "cm"), colour = "black"),axis.ticks = element_line(colour = "black")) +
    scale_colour_manual(name = "Treatment",
                        labels = trt_name,
                        values=c("#028ff7","#f76a02", "#8ecbf8", "#fec398")) +   
    scale_shape_manual(name = "Treatment",
                       labels = trt_name,
                       values = c(19, 19, 17, 17))
  print(p)
}

Fig 2a Bacteria initial growth

trt_name <- c(expression(paste("Full 22 ",degree,"C")), expression(paste("Full 25 ",degree,"C")), expression(paste("Half 22 ",degree,"C")),expression(paste("Half 25 ",degree,"C")) )

# Bacteria initial growth
TEB_interactive(bac_desc,bac_desc$nut,bac_desc$init_growth,bac.mean_desc,bac.mean_desc$nut,bac.mean_desc$ig,"Initial Growth")

Fig 2b Tetrahymena initial growth

# Tetrahymena initial growth
TEB_interactive(tetra_desc,tetra_desc$nut,tetra_desc$tetra.ig,t.mean_desc,t.mean_desc$nut,t.mean_desc$ig,"Initial Growth")

Fig 2c Euplotes initial growth

# Euplotes initial growth
TEB_interactive(eup_desc,eup_desc$nut,eup_desc$init_growth,eup.mean_desc,eup.mean_desc$nut,eup.mean_desc$ig,"Initial Growth")

Fig 2 d-f Maximum abundance

Fig 2d Bacteria Maximum Abundance

TEB_interactive2 <- function(data1,x1,y1,data2,x2,y2,y_label){
  p <- ggplot(data=data1, aes(x=x1, y=y1, group=temp, shape= tre,colour=tre)) +
     geom_line(data = data2, aes(x=x2, y=y2), size = 0.6, color = "#c3c3c3") +
    scale_x_discrete(limits=c("Half","Full")) +
    labs(y=y_label, x = "Nutrient") +
    geom_jitter(size = 2, shape = rep(c(1,2),c(36,36 )), width = 0.25) +
    stat_summary(fun.data = mean_se, geom = "errorbar", width=0.2, size = 0.5) +
    geom_point(data = data2, aes(x=x2, y=y2, group=temp, shape=tre, colour=tre), size = 3) +
    theme_bw() +
   theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),axis.ticks.length=unit(-0.12, "cm"), axis.text.x = element_text(margin=unit(c(0.25,0.25,0.25,0.25), "cm"),colour = "black"), axis.text.y = element_text(margin=unit(c(0.25,0.25,0.25,0.25), "cm"), colour = "black"),axis.ticks = element_line(colour = "black")) +
    scale_colour_manual(name = "Treatment",
                        labels = trt_name,
                        values=c("#028ff7","#f76a02", "#8ecbf8", "#fec398")) +   
    scale_shape_manual(name = "Treatment",
                       labels = trt_name,
                       values = c(19, 19, 17, 17))
  print(p)
}

# Bacteria Maximum Abundance
TEB_interactive2(bac_max,bac_max$nut,bac_max$OD,bac_max_mean,bac_max_mean$nut,bac_max_mean$ODmax,"Maximum Abundance")

Fig 2e Tetrahymena Maximum Abundance

# Tetrahymena Maximum Abundance
TEB_interactive(tetra_desc,tetra_desc$nut,tetra_desc$max,t.mean_desc,t.mean_desc$nut,t.mean_desc$max.m,"Maximum Abundance")

Fig 2f Euplotes Maximum Abundance

# Euplotes Maximum Abundance
TEB_interactive(eup_desc,eup_desc$nut,eup_desc$max,eup.mean_desc,eup.mean_desc$nut,eup.mean_desc$max.m,"Maximum Abundance")

Fig 2 g-i Coefficient of vairation (CV)

Fig 2g Bacteria CV

# Bacteria CV
TEB_interactive(bac_desc,bac_desc$nut,bac_desc$CV,bac.mean_desc,bac.mean_desc$nut,bac.mean_desc$CV.m,"CV")

Fig 2h Tetrahymena CV

# Tetrahymena CV
TEB_interactive(tetra_desc,tetra_desc$nut,tetra_desc$CV,t.mean_desc,t.mean_desc$nut,t.mean_desc$CV.m,"CV")

Fig 2i Euplotes CV

# Euplotes CV
TEB_interactive(eup_desc,eup_desc$nut,eup_desc$CV,eup.mean_desc,eup.mean_desc$nut,eup.mean_desc$CV.m,"CV")

Fig 2j Day of collapse of Tetrahymena

# Tetrahymena Time Collapes
TEB_interactive(tetra_desc,tetra_desc$nut,tetra_desc$tc,t.mean_desc,t.mean_desc$nut,t.mean_desc$tc.m,"Day Collapsed") 

Fig 2k Day of peak of Euplotes

 # Euplotes Day peaking
TEB_interactive(eup_desc,eup_desc$nut,eup_desc$daypeak,eup.mean_desc,eup.mean_desc$nut,eup.mean_desc$daypeak.m,"Day Peaked")

Fig 3 Significant effects of Temperature and nutrients

All the linear models for temperature and nutrients effects can be found in appendix I Table S3 and the code for the models can be found in Rmarkdown file for appendix. Here we added up the significant additive temperature and nutrients effects and their interactions on species initial growh, maximum abundance and CV to plot Figure 3.

TN_count <- as.data.frame(matrix(c(rep(NA,12) ), 3, 4))
colnames(TN_count) <- c("species","Temp", "Nut","T*N")
TN_count[,1] <- c("Bacteria", "Tetrahymena","Euplotes")
TN_count[1,2:4] <- c(1, 2, 0)
TN_count[2,2:4] <- c(3, 2, 2)
TN_count[3,2:4] <- c(2, 3, 2)
TN_count_ga <-gather(TN_count,key=effects, value = Counts,"Temp", "Nut","T*N" );TN_count_ga
 ggplot(data=TN_count_ga, aes(x=species, y=Counts, fill=factor(effects,levels = c("T*N","Nut","Temp")))) +
  geom_bar(stat="identity", width=0.7)+ 
  scale_x_discrete(limits = c("Bacteria", "Tetrahymena","Euplotes"))+
   scale_y_discrete(limits = c(0,3,6,9))+
  scale_fill_manual(name = "Effects",labels = c("N*T","N", "T"),values=c("#efd3d7" ,"#cbc0d3","#8e9aaf")) +   
  theme_bw()+theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), axis.ticks.length=unit(-0.2, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")))

Fig 4 CCM predictability plot between speceices ecological and phenotypic dynamcis

4.1 Interpolation Using Three Different Methods

4.1.1 approx (linear interpolation)

4.1.1.1 Density data

Here we generate linear interpolation prediction for the density of all three species. First, we defined a function to interpolate by treatment and by replication.

linearapprox <- function(data,Yvar, tre.vector, rep.vector,daycolnum,ycolnum,x_out){
  total_ap <- data.frame()
  for (i in tre.vector){
  for(j in rep.vector){
  ap <- as.data.frame(approx(data$day[which(data$treatNum== i & data$rep == j)], Yvar[which(data$treatNum== i & data$rep == j)], xout = x_out,method = "linear", rule = 1, f= 0))
d <- data[which(data$treatNum== i & data$rep == j),c(daycolnum,ycolnum)] 
colnames(d) <- c("x", "y")
t<- rbind(ap,d) 
t$treatNum <- rep(c(i), n = 16)
t$rep <- rep(c(j), n = 16)
t <- t %>% arrange(x)
mCCM_nas <- as.data.frame(matrix(rep(c(16,NA,i,j,17,NA,i,j)),byrow = T,nrow = 2 )) # NAs added here for multispcaial CCM analysis formatting. Each time series should be spererated by two NAs. See instruction for "SSR_pred_boot" function
colnames(mCCM_nas) <- c("x", "y", "treatNum","rep")
t2 <- rbind(t,mCCM_nas)
total_ap <- rbind(total_ap, t2)
    }
  }
colnames(total_ap) <- c("day", "density", "treatNum","rep")
total_ap$density[total_ap$density < 0] = 0  #***** all the negative prediction are replace by 0!
print(total_ap)
}

We then interpolate all the missing data.

tetra_approx <- linearapprox(tetra, tetra$density, c(1:4), c(1:6), 3,8, c(5:6, 12:13)) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep),species = c(rep("tetra", 432)))
eup_approx <- linearapprox(eup, eup$density, c(1:4), c(1:6), 3,8,c(5:6, 12:13)) %>% mutate(treatNum = as.factor(treatNum),  rep = as.factor(rep),species = c(rep("eup", 432)))
bac_approx <- linearapprox(bac, bac$OD, c(1:4), c(1:5), 4,9,c(5:6, 12:13)) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep),species = c(rep("bac", 339)))
# all species
all_approx <- rbind(tetra_approx, eup_approx, bac_approx)
# plot
trt_name = c(expression(paste("Full22 ",degree,"C")),expression(paste("Full25 ",degree,"C")), expression(paste("Half22 ",degree,"C")),expression(paste("Half25 ",degree,"C")))

Here are the ecological time serieses with linear interpolation for all species.

ggplot(data = all_approx, aes(x=day, y=density, group=treatNum, color=treatNum, shape= treatNum)) +
    geom_point(size =3) +
facet_wrap(~species, scale = "free")+
  scale_colour_manual(name = "Treatment",
                        labels = trt_name,
                        values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment",
                       labels = trt_name,
                       values = c(19, 19, 17, 17))

##### 4.1.1.2 Trait data

Now we use linear interpolation to predict the trait data of both protist species on missing days. Because of the low density of Euplotes in the beginning of the experiment, there are days where we detected no Euplotes in our samples, thus no trait data were recorded. Here we also interpolate trait values for those days to prevent NAs in the time series (multispacial CCM) analysis later on.

# Interpolation for the missing days, day 5,6,12,13
t.t_app <- linearapprox(tetraSumm, tetraSumm$area, c(1:4), c(1:6), 4,15,c(5,6,12,13)) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep),species = c(rep("tetra", 432)))
colnames(t.t_app) <- c("day","area", "treatNum", "rep","species" )

e.t_app <- linearapprox(eupSumm, eupSumm$area, c(1:4), c(1:6), 4,15,c(5,6,12,13)) %>% mutate(treatNum = as.factor(treatNum),  rep = as.factor(rep),species = c(rep("eup", 432)))
colnames(e.t_app) <- c("day","area", "treatNum", "rep","species" )

#  Selecting missing trait data of Euplotes for linear interpolation
e.t.na <- e.t_app[is.na(e.t_app$area),c(1,3,4)] %>% filter(!day %in% c(17,18)) # row for day 17, 18 were added for multispacial CCM analysis data format only (each time series is seperated by two NAs).Details see instruction for "SSR_pred_boot" function

# interpolation missing trait data of Euplotes
for (i in 1:nrow(e.t.na)){
  d <- approx(eupSumm$day[which(eupSumm$treatNum== e.t.na[i,2] & eupSumm$rep == e.t.na[i,3])], eupSumm$area[which(eupSumm$treatNum== e.t.na[i,2] & eupSumm$rep == e.t.na[i,3])], xout = e.t.na[i,1], method = "linear", rule = 1, f= 0)
  e.t_app$area[which(e.t_app$treatNum == e.t.na[i,2] & e.t_app$rep == e.t.na[i,3] & e.t_app$day == e.t.na[i,1])] <- d$y
}

Here are T.pyriformis and Eupltoes sp. phenotypic dynamics with linear interpolation.

trait_app <- rbind(t.t_app, e.t_app)
ggplot(data = trait_app, aes(x=day, y=area, group=treatNum, color=treatNum, shape= treatNum)) +
    geom_point(size =3) +
facet_wrap(~species, scale = "free")+
  scale_colour_manual(name = "Treatment",
                        labels = trt_name,
                        values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment",
                       labels = trt_name,
                       values = c(19, 19, 17, 17))

4.1.2 spline

4.1.2.1 Density data

As before, here we define function to interpolate missing data with spline.

require(stats)
intpl_spline <- function(data,Yvar, tre.vector, rep.vector,daycolnum,denscolnum){
  total_sp <- data.frame()
  for (i in tre.vector){
  for(j in rep.vector){
  sp1 <- as.data.frame(spline(data$day[which(data$treatNum== i & data$rep == j)], Yvar[which(data$treatNum== i & data$rep == j)], method ="fmm", n = 2,xmin = 5, xmax = 6))
  sp2 <- as.data.frame(spline(data$day[which(data$treatNum== i & data$rep == j)], Yvar[which(data$treatNum== i & data$rep == j)], method ="fmm", n = 2,xmin = 12, xmax = 13))
d <- data[which(data$treatNum== i & data$rep == j),c(daycolnum,denscolnum)] 
colnames(d) <- c("x", "y")
t<- rbind(sp1,d) %>% rbind(.,sp2)
t$treatNum <- rep(c(i), n = 16)
t$rep <- rep(c(j), n = 16)
t <- t %>% arrange(x)
mCCM_nas <- as.data.frame(matrix(rep(c(16,NA,i,j,17,NA,i,j)),byrow = T,nrow = 2 ))
colnames(mCCM_nas) <- c("x", "y", "treatNum","rep")
t2 <- rbind(t,mCCM_nas)
total_sp <- rbind(total_sp, t2)
    }
  }
colnames(total_sp) <- c("day", "density", "treatNum","rep")
total_sp$density[total_sp$density < 0] = 0  #***** all the negative prediction are replace by 0!
print(total_sp)
}

Then we predict all the missing data in species ecological dynamics as in section 4.1.1 but with spline interpolation.

tetra_spline <- intpl_spline(tetra, tetra$density, c(1:4), c(1:6), 3,8) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep),species = c(rep("tetra", 432)))
eup_spline <- intpl_spline(eup, eup$density, c(1:4), c(1:6), 3,8) %>% mutate(treatNum = as.factor(treatNum),  rep = as.factor(rep),species = c(rep("eup", 432)))
bac_spline <- intpl_spline(bac, bac$OD, c(1:4), c(1:5), 4,9) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep),species = c(rep("bac", 339)))
# all species
all_spline <- rbind(tetra_spline, eup_spline, bac_spline)

Here are the species ecological dynamics with spline interpolation.

ggplot(data = all_spline, aes(x=day, y=density, group=treatNum, color=treatNum, shape= treatNum)) +
    geom_point(size =3) +
facet_wrap(~species, scale = "free")+
  scale_colour_manual(name = "Treatment",
                        labels = trt_name,
                        values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment",
                       labels = trt_name,
                       values = c(19, 19, 17, 17))

4.1.2.2 Trait data

As in section 4.1.1.2, we interpolate the missing trait data of T. pyriformis and Euplotes sp. but with spline interpolation.

# Interpolated the trait data on the missing days with previously defined function.
t.t_sp <- intpl_spline(tetraSumm, tetraSumm$area, c(1:4), c(1:6), 4,15) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep),species = c(rep("tetra", 432)))
colnames(t.t_sp) <- c("day","area", "treatNum", "rep","species" )
tetraSumm[which(tetraSumm$day == c(11,14)), c(2,4,15)]
t.t_sp$area[which(t.t_sp$day== c(12,13))] <- NA # because all treatments had average extinction time for tetrahymena before or day 11, that means tetrahymena is functionally extincted within the microcosm. Spline() predicted trait/area value on day 12-13 is not reliable as they are negative  while there is only one microcosm out of 24 had trait measure on day 11 and none on day 14. There for I replaced them with NAs.


e.t_sp <- intpl_spline(eupSumm, eupSumm$area, c(1:4), c(1:6), 4,15) %>% mutate(treatNum = as.factor(treatNum),  rep = as.factor(rep),species = c(rep("eup", 432)))
colnames(e.t_sp) <- c("day","area", "treatNum", "rep","species" )
e.t_sp$area

# Interpolate missing trait data of Euplotes.
e.t.na <- e.t_sp[is.na(e.t_sp$area),c(1,3,4)] %>% filter(!day %in% c(16,17)) # rows for day 16,17 were added for multispacial CCM analysis data format only (each time series is seperated by two NAs). 

for (i in 1:nrow(e.t.na)){
  d <- spline(eupSumm$day[which(eupSumm$treatNum== e.t.na[i,2] & eupSumm$rep == e.t.na[i,3])], eupSumm$area[which(eupSumm$treatNum== e.t.na[i,2] & eupSumm$rep == e.t.na[i,3])], method ="fmm", n = 2,xmin = e.t.na[i,1], xmax = e.t.na[i,1])
  e.t_sp$area[which(e.t_sp$treatNum == e.t.na[i,2] & e.t_sp$rep == e.t.na[i,3] & e.t_sp$day == e.t.na[i,1])] <- d$y[1]
}
e.t_sp$area

Here are T.pyriformis and Eupltoes sp. phenotypic dynamics with Spline interpolation.

trait_sp <- rbind(t.t_sp, e.t_sp)
ggplot(data = trait_sp, aes(x=day, y=area, group=treatNum, color=treatNum, shape= treatNum)) +
    geom_point(size =3) +
facet_wrap(~species, scale = "free")+
  scale_colour_manual(name = "Treatment",
                        labels = trt_name,
                        values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment",
                       labels = trt_name,
                       values = c(19, 19, 17, 17))

4.1.3 smooth.spline

4.1.3.1 Density data

Again, we first define function that allow us to interpolate missing data by treatment and by repliate with smooth spline interpolation.

intpl_smooth_spline <- function(data,Yvar, tre.vector, rep.vector,daycolnum,denscolnum,x_out){
  total_sp <- data.frame()
  for (i in tre.vector){
  for(j in rep.vector){
  sp1 <- smooth.spline(data$day[which(data$treatNum== i & data$rep == j)], Yvar[which(data$treatNum== i & data$rep == j)])
  p <- predict(sp1,x_out)
d <- data[which(data$treatNum== i & data$rep == j),c(daycolnum,denscolnum)] 
colnames(d) <- c("x", "y")
t<- rbind(p,d)
t$treatNum <- rep(c(i), n = 16)
t$rep <- rep(c(j), n = 16)
t <- t %>% arrange(x)
mCCM_nas <- as.data.frame(matrix(rep(c(16,NA,i,j,17,NA,i,j)),byrow = T,nrow = 2 ))
colnames(mCCM_nas) <- c("x", "y", "treatNum","rep")
t2 <- rbind(t,mCCM_nas)
total_sp <- rbind(total_sp, t2)
    }
  }
colnames(total_sp) <- c("day", "density", "treatNum","rep")
total_sp$density[total_sp$density < 0] = 0  #***** all the negative prediction are replace by 0!
print(total_sp)
}

Now we interpolation missing of species ecological time series.

tetra_sm_sp <- intpl_smooth_spline(tetra, tetra$density, c(1:4), c(1:6), 3,8,c(5,6,12,13)) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep),species = c(rep("tetra", 432))) 
eup_sm_sp <- intpl_smooth_spline(eup, eup$density, c(1:4), c(1:6), 3,8,c(5,6,12,13)) %>% mutate(treatNum = as.factor(treatNum),  rep = as.factor(rep),species = c(rep("eup", 432)))
bac_sm_sp <- intpl_smooth_spline(bac, bac$OD, c(1:4), c(1:5), 4,9,c(5,6,12,13)) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep),species = c(rep("bac", 339)))
all_sm_sp <- rbind(tetra_sm_sp,eup_sm_sp, bac_sm_sp)

Here are the ecological time series of all species with smooth spline interpolation.

ggplot(data = all_sm_sp, aes(x=day, y=density, group=treatNum, color=treatNum, shape= treatNum)) +
    geom_point(size =3) +
facet_wrap(~species, scale = "free")+
  scale_colour_manual(name = "Treatment",
                        labels = trt_name,
                        values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment",
                       labels = trt_name,
                       values = c(19, 19, 17, 17))

4.1.3.2 Trait data

As in section 4.1.1.2, we interpolate the missing data but with smooth spline interpolation.

tetraSumm.no.na <- tetraSumm %>% na.omit()
t.t_smsp <- intpl_smooth_spline(tetraSumm.no.na, tetraSumm.no.na$area, c(1:4), c(1:6), 4,15, c(5,6)) %>% mutate(treatNum = as.factor(treatNum), rep = as.factor(rep))
t.t_smsp$species = rep(c("tetra"), times = nrow(t.t_smsp))
colnames(t.t_smsp) <- c("day","area", "treatNum", "rep","species" )
t.t_smsp$area

# only interpolated the trait for the missing days
eupSumm.no.na <- eupSumm %>% na.omit()
e.t_smsp0 <- intpl_smooth_spline(eupSumm.no.na, eupSumm.no.na$area, c(1:4), c(1:6), 4,15, c(5,6,12,13)) %>% mutate(treatNum = as.factor(treatNum),  rep = as.factor(rep),species = c(rep("eup", nrow(.))))
colnames(e.t_smsp0) <- c("day","area", "treatNum", "rep","species" )
e.t_smsp <- e.t.na %>% mutate(area = rep(NA, times= nrow(e.t.na)), species = rep("eup", times= nrow(e.t.na))) %>% rbind(.,e.t_smsp0) %>% arrange(treatNum, rep, day)
e.t_smsp$area 

# Because of the low density of Euplotes in the beginning of the experiment, there are days where we detected no Euplotes in our samples, thus no trait data. 
# I also interpolated trait values for days that we didn't detected individual for trait measurements. 
e.t.na.smsp <- e.t.na %>% mutate(area = rep(NA, times= nrow(e.t.na)), species = rep("eup", times= nrow(e.t.na)))# missing values in area
for (i in 1:nrow(e.t.na.smsp)){
  sp <- smooth.spline(eupSumm.no.na$day[which(eupSumm.no.na$treatNum== e.t.na[i,2] & eupSumm.no.na$rep == e.t.na[i,3])], eupSumm.no.na$area[which(eupSumm.no.na$treatNum== e.t.na[i,2] & eupSumm.no.na$rep == e.t.na[i,3])])
  d <- predict( sp, e.t.na[i,1])
  e.t.na.smsp$area[i] <- d$y
}

e.t_smsp <- rbind(e.t_smsp0,e.t.na.smsp) %>% arrange(treatNum, rep,day)
e.t_smsp$area

Here are T.pyriformis and Eupltoes sp. phenotypic dynamics with Smooth Spline interpolation.

trait_smsp <- rbind(t.t_smsp, e.t_smsp)
ggplot(data = trait_smsp, aes(x=day, y=area, group=treatNum, color=treatNum, shape= treatNum)) +
    geom_point(size =3) +
facet_wrap(~species, scale = "free")+
  scale_colour_manual(name = "Treatment",
                        labels = trt_name,
                        values=c("#028ff7","#f76a02", "#8ecbf8", "#fdc470")) +   
    scale_shape_manual(name = "Treatment",
                       labels = trt_name,
                       values = c(19, 19, 17, 17))

4.2 Multispacial CCM

4.2.1 Using linear interpolation for Multispacial CCM

Here are the ecological data with linear interpolation that we will use for multispacial CCM. We first spread the ecological data so it’s easier to index for multispacial CCM analysis.

tetra_app_ccm <- tetra_approx %>%spread( key = treatNum, value = density) %>% arrange(rep,day) %>% rename(F22 = 4, F25 = 5, H22 = 6, H25 =7)
eup_app_ccm  <- eup_approx %>%spread( key = treatNum, value = density) %>% arrange(rep,day)%>% rename(F22 = 4, F25 = 5, H22 = 6, H25 =7)
bac_app_ccm  <- bac_approx %>%spread( key = treatNum, value = density) %>% arrange(rep,day)%>% rename(F22 = 4, F25 = 5, H22 = 6, H25 =7)

# remove day 15 for tetra and eup to run multispacial CCM with bacteria since Day 15 data is missing for bacteria
    tetra_app_ccm_d14 <- tetra_app_ccm[-which(tetra_app_ccm$day == 15),]
    eup_app_ccm_d14 <- eup_app_ccm[-which(tetra_app_ccm$day == 15),]

Here we calculated the optimal embeding dimension for each of the species by treatment. First, we define function to calculate and record optimal embedding dimension by treatment for multispacial CCM analysis.

opt_E <- function(maxE, data, Name){
  Emat<-as.data.frame(matrix(nrow=maxE-1, ncol=5)); colnames(Emat)<-c("Embed_D", "F22", "F25", "H22","H25"); Emat$Embed_D <- c(2:maxE)
  for(E in 2:maxE) {
    for (j in 1:4){
  Emat[E-1,j+1]<-SSR_pred_boot(A=data[,j+3], E=E, predstep=1, tau=1)$rho
  }}
 optE <- Emat %>% pivot_longer(!Embed_D,names_to = "trt", values_to = "rho" )%>% group_by(trt) %>% filter(rho == max(rho)) %>% select(trt,Embed_D) %>% as.data.frame()
 colnames(optE) <- c("trt", Name)
 return(optE)
}

Then we calculate the optimal embedding dimension for each species by treatment.

bac_app_E <- opt_E(6,bac_app_ccm, "bac")
tetra_app_E <- opt_E(6,tetra_app_ccm, "tetra")
eup_app_E <- opt_E(6,eup_app_ccm, "eup")
tetra_app_E_d14 <- opt_E(6,tetra_app_ccm_d14, "tetra_d14")
eup_app_E_d14 <- opt_E(6,eup_app_ccm_d14, "eup_d14")
opt_E_app1 <- left_join(bac_app_E,tetra_app_E, by = "trt") %>% left_join(eup_app_E, by = "trt") %>% left_join(tetra_app_E_d14, by = "trt") %>% left_join(eup_app_E_d14, by = "trt") %>% arrange(trt)

Here are species trait data with linear interpolation prediction that will be used in multispacial CCM.

t.t_app_ccm <- t.t_app %>%spread( key = treatNum, value = area) %>% arrange(rep,day) 
    colnames(t.t_app_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")

t.t_app_ccm_d9 <- t.t_app_ccm %>% filter(! day %in% c(10:15 ))# only kept tetrahymena trait data from day 0 - day 9 to avoid problems with too many NAs in the time series for ccm test. On average tetrahymena population went extinct by day 9.

e.t_app_ccm <- e.t_app %>%spread( key = treatNum, value = area) %>% arrange(rep,day) # only kept tetrahymena trait data from day 0 - day 9 since on average tetrahymena population went extinct by day 9.
    colnames(e.t_app_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")
e.t_app_ccm_d9 <- e.t_app_ccm %>% filter(! day %in% c(10:15 ))

# crop ecological dynamics that are interpolated by linear interpolation (approx)
bac_app_ccm_d9 <- bac_app_ccm %>% filter(! day %in% c(10:14))
tetra_app_ccm_d9 <- tetra_app_ccm %>% filter(! day %in% c(10:14))
eup_app_ccm_d9 <- eup_app_ccm %>% filter(! day %in% c(10:14))

Calculate the best embeding dimension for the trait time series, and any cropped time series to match the length of tetrahymena trait time series. Time series’s embeding dimension that has not been calculated are: - 1. tetrahymena trait (cropped to day 9) - 2. Euplotes trait - 3. Euplotes trait crop to day 9 - 4-6. bacteria, tetrahymena, euplotes ecological dynamics/ density time series crop to day 9

t.t_app_E <- opt_E(6,t.t_app_ccm_d9, "t.t")
e.t_app_E <- opt_E(6,e.t_app_ccm, "e.t")
e.t.d9_app_E <- opt_E(6,e.t_app_ccm_d9, "e.t_d9")
bac.d9_app_E <- opt_E(6,bac_app_ccm_d9, "bac_d9")
tetra.d9_app_E <- opt_E(6,tetra_app_ccm_d9, "tetra_d9")
eup.d9_app_E <- opt_E(6,eup_app_ccm_d9, "eup_d9")
opt_E_app <- left_join(opt_E_app1, t.t_app_E) %>% left_join(e.t_app_E, by = "trt") %>% left_join(e.t.d9_app_E, by = "trt") %>% left_join(bac.d9_app_E, by = "trt") %>% left_join(tetra.d9_app_E, by = "trt") %>% left_join(eup.d9_app_E, by = "trt")  %>% arrange(trt)

Here is a list of all the dataframes that will be used for multispacial CCM test.

dlist_app <- list(bac_app_ccm, tetra_app_ccm,tetra_app_ccm_d14,eup_app_ccm,eup_app_ccm_d14, t.t_app_ccm, e.t_app_ccm,e.t_app_ccm_d9,bac_app_ccm_d9, tetra_app_ccm_d9,eup_app_ccm_d9)
names(dlist_app) <- c("bac_app_ccm", "tetra_app_ccm", "tetra_app_ccm_d14", "eup_app_ccm", "eup_app_ccm_d14", "t.t_app_ccm", "e.t_app_ccm","e.t_app_ccm_d9", "bac_app_ccm_d9", "tetra_app_ccm_d9", "eup_app_ccm_d9")
4.2.1.1 Multispacial CCM test

Here we define the function that runs multispcial CCM test on designated pairs of times series.

msCCM_list <- function(dlist,testinfo, optE, itrtn_num){
msCCM_data <- NULL
trename <- c("F22", "F25", "H22", "H25")
for(j in 1:4){
  ccmlist <- NULL
  for(i in 1: nrow(testinfo)) {
  ccmlist[[testinfo[i,"pair"]]] <- CCM_boot(dlist[[testinfo[i,"cause"]]][,j+3], dlist[[testinfo[i,"reciever"]]][,j+3], optE[j,testinfo[i,"E"]], tau=1, iterations = itrtn_num)
  }
  lname <- names(ccmlist)
  names(ccmlist) <- paste(rep(trename[j],nrow(testinfo)),lname, sep = "_")
  msCCM_data<- c(msCCM_data,ccmlist)
}
return(msCCM_data)
}

Eco-Eco interactions. Here we performed the multispacial CCM for the ecological dynamics time series for each of the species pair by treatment. Because we are testing the causal relationship between each pair, we ran the multispacial CCM test ‘CCMboot()’ on each species pair twice, with different order in the CCMboot() test. This results 2 casual relationship test per species pair * 3 species pair * 4 treatments = 24 tests results. Each of the test run 800 iteration.

test_app_EcEc <- NULL # creating a dataframe to call timeseries for the CCM tests, loop for each treatments
test_app_EcEc$pair <- c("b.t", "t.b","b.e","e.b","t.e","e.t")
test_app_EcEc$cause <- c("bac_app_ccm","tetra_app_ccm_d14", "bac_app_ccm","eup_app_ccm_d14","tetra_app_ccm","eup_app_ccm")
test_app_EcEc$reciever <-  c("tetra_app_ccm_d14","bac_app_ccm", "eup_app_ccm_d14","bac_app_ccm","eup_app_ccm","tetra_app_ccm")
test_app_EcEc$E <- c("bac","tetra_d14", "bac","eup_d14","tetra","eup")
test_app_EcEc <- test_app_EcEc %>% as.data.frame()
msCCM_approx_EcEc <- msCCM_list(dlist_app, test_app_EcEc, opt_E_app, 800)

Eco-Pheno interactions. Here we use the multispacial CCM test to examine the causal effects of the bacterial and protists ecological dynamics time series on protists phenotypic dynamics time series (the changes in protists size over time).

test_app_EcPhn <- NULL # creating a dataframe to call time series for the CCM tests, loop for each treatments
test_app_EcPhn$pair <- c("b.tt", "b.et","t.tt","t.et","e.tt","e.et")
test_app_EcPhn$cause <- c("bac_app_ccm_d9","bac_app_ccm", "tetra_app_ccm_d9","tetra_app_ccm","eup_app_ccm_d9","eup_app_ccm")
test_app_EcPhn$reciever  <- c("t.t_app_ccm","e.t_app_ccm", "t.t_app_ccm","e.t_app_ccm","t.t_app_ccm","e.t_app_ccm")
test_app_EcPhn$E <- c("bac_d9","bac", "tetra_d9","tetra","eup_d9","eup")
test_app_EcPhn <- test_app_EcPhn %>% as.data.frame()               
msCCM_approx_EcPhn <- msCCM_list(dlist_app, test_app_EcPhn, opt_E_app, 800)

Pheno-Eco interactions. Now we examine the causal effects of of the size of both protists species on the ecological dynamcis of all three species.

test_app_PhnEc <- NULL # creating a dataframe to call time series for the CCM tests, loop for each treatments
test_app_PhnEc$pair <- c("tt.b", "et.b","tt.t","et.t","tt.e","et.e")
test_app_PhnEc$cause <- c("t.t_app_ccm","e.t_app_ccm", "t.t_app_ccm","e.t_app_ccm","t.t_app_ccm","e.t_app_ccm")
test_app_PhnEc$reciever  <- c("bac_app_ccm_d9","bac_app_ccm", "tetra_app_ccm_d9","tetra_app_ccm","eup_app_ccm_d9","eup_app_ccm")
test_app_PhnEc$E <- c("t.t", "e.t", "t.t", "e.t", "t.t", "e.t")
test_app_PhnEc <- test_app_PhnEc %>% as.data.frame()                          
msCCM_approx_PhnEc <- msCCM_list(dlist_app, test_app_PhnEc, opt_E_app, 800)

Pheno-Pheno interactions.

test_app_PhnPhn <- NULL 
test_app_PhnPhn$pair <- c("tt.et", "et.tt")
test_app_PhnPhn$cause <- c("t.t_app_ccm","e.t_app_ccm_d9")
test_app_PhnPhn$reciever  <- c("e.t_app_ccm_d9", "t.t_app_ccm")
test_app_PhnPhn$E <- c("t.t", "e.t_d9")
test_app_PhnPhn <- test_app_PhnPhn %>% as.data.frame()                          
msCCM_approx_PhnPhn <- msCCM_list(dlist_app, test_app_PhnPhn, opt_E_app, 800)
4.2.1.2 Multispactial CCM plots with Linear interpolation (Apendix II Fig 1-4)

Here we define a new function to plot multispacial CCM predictability against library length.

plot_msccm <- function (data, i, Title, color, linestyle){
  plot(data[[i]][["Lobs"]], data[[i]][["rho"]],type="l", col = color[1], lwd = 2, lty = linestyle[1], xlim=c(5, 100), ylim=c(0,1), xlab="Library length", ylab="rho",cex.lab=1.5, cex.axis=1.5, cex.names = 1.5 ) # b.t
  lines(data[[i+1]][["Lobs"]], data[[i+1]][["rho"]], type="l", col = color[2], lwd = 2, lty = linestyle[2]) 
  lines(data[[i+2]][["Lobs"]], data[[i+2]][["rho"]], type="l", col = color[3], lwd = 2, lty = linestyle[3]) 
  lines(data[[i+3]][["Lobs"]], data[[i+3]][["rho"]], type="l", col = color[4], lwd = 2, lty = linestyle[4])
  lines(data[[i+4]][["Lobs"]], data[[i+4]][["rho"]],type="l", col= color[5], lty=linestyle[5],lwd=1.5) 
  lines(data[[i+5]][["Lobs"]], data[[i+5]][["rho"]],col = color[6], lty=linestyle[6], lwd=1.5)
  title(main = Title, cex.main = 2)}
Appendix II Fig 1
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
# Eco-Eco msCCM
EcEc_col <- c( "#A997DF", "#A997DF","#add0af","#add0af","#FFC857","#FFC857")
EcEc_line <- c(1,2,1,2,1,2)
{for(i in c(1:4)){
    plot_msccm(msCCM_approx_EcEc, i*6-5,trt_name[i], EcEc_col,EcEc_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = EcEc_col, lty = EcEc_line, legend =c("B -> T","T -> B", "B -> E", "E -> B", "T -> E","E -> T" ), lwd = 2,  cex = 2)
mtext("Eco-Eco multispacial CCM Prediction (linear)", outer = TRUE, cex = 2)}

Appendix II Fig 2
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
# Eco-Pheno msCCM 
EcPhn_col <- c("#A997DF","#add0af","#929084","#FFC857","#FFC857","#2E4052")
EcPhn_line <- c(1,1,3,1,2,3)
{for(i in c(1:4)){
    plot_msccm(msCCM_approx_EcPhn, i*6-5,trt_name[i],EcPhn_col,EcPhn_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = EcPhn_col, lty = EcPhn_line, legend =c( "B -> Tt", "B -> Et", "T -> Tt","T -> Et", "E -> Tt", "E -> Et"), lwd = 2,  cex = 2)
mtext("Eco-Pheno multispacial CCM Prediction (linear)", outer = TRUE, cex = 2)}

Appendix II Fig 3
# Pheno-Eco msCCM 
PhnEc_col <- c("#A997DF","#add0af","#929084","#FFC857","#FFC857","#2E4052")
PhnEc_line <- c(2,2,3,2,1,3)

par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
{for(i in c(1:4)){
    plot_msccm(msCCM_approx_PhnEc, i*6-5,trt_name[i],PhnEc_col, PhnEc_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = PhnEc_col, lty = PhnEc_line, legend =c( "Tt -> B", "Et -> B", "Tt -> T","Et -> T", "Tt -> E", "Et -> E"), lwd = 2,  cex = 2)
mtext("Pheno-Eco multispacial CCM Prediction (linear)", outer = TRUE, cex = 2)}

Appendix II Fig 4
#Pheno-Pheno
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
for(i in c(1,3,5,7)){
plot(msCCM_approx_PhnPhn[[i]][["Lobs"]], msCCM_approx_PhnPhn[[i]][["rho"]],type="l", col = "#FFC857", lwd = 2, lty = 1, xlim=c(5, 80), ylim=c(-0,1), xlab="Library length", ylab="rho",cex.lab=1.5, cex.axis=1.5, cex.names = 1.5 ) # b.t
  lines(msCCM_approx_PhnPhn[[i+1]][["Lobs"]], msCCM_approx_PhnPhn[[i+1]][["rho"]], type="l", col = "#FFC857", lwd = 2, lty = 2) 
}
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = c("#FFC857","#FFC857"), lty = c(1,2), legend =c( "Tt -> Et", "Et -> Tt"), lwd = 2,  cex = 2)
mtext("Pheno-Pheno multispacial CCM Prediction (linear)", outer = TRUE, cex = 2)

4.2.1.3 Plot the rho at largest library size

Here I set up the data frame to store the multispacial CCM results, and then I extract the predictability (rho) at the largest library length of the each multispcial CCM test.

# Eco-Eco
EcEcRho<- str_split_fixed(names(msCCM_approx_EcEc), '_', 2) %>% as.data.frame %>% rename(trt = V1,mapping = V2) %>% mutate(temp = rep(c("22C", "25C"), each = 6, time = 2),nut=rep(c("High", "Low"), each = 12), species = rep(c("BT", "BE", "TE"),each=2, time = 4), direc = rep(c("up", "down"), time = 12), dynamics = rep(c("EcoEco"), time = 24))

Extra.rho <- function(i,dlist,factor) {last(dlist[[i]]$rho)}
Extra.Lobs <- function(i,dlist,factor) {last(dlist[[i]]$Lobs)}

EcEcRho$rho_app <- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_approx_EcEc))
EcEcRho$Lobs_app <- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_approx_EcEc))

# Eco-Pheno
EcPhnRho<- str_split_fixed(names(msCCM_approx_EcPhn), '_', 2) %>% as.data.frame %>% rename(trt = V1,mapping = V2) %>% mutate(temp = rep(c("22C", "25C"), each = 6, time = 2),nut=rep(c("High", "Low"), each = 12),species = rep(c("BT", "BE", "Tetra","TE","TE","Eup"),time = 4), direc = rep(c("up", "up","single","up","down", "single"), time = 4), dynamics = rep(c("EcoPheno"), time = 24))
EcPhnRho$rho_app <- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_approx_EcPhn ))
EcPhnRho$Lobs_app <- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_approx_EcPhn))

# Pheno-Eco
PhnEcRho<- str_split_fixed(names(msCCM_approx_PhnEc), '_', 2) %>% as.data.frame %>% rename(trt = V1,mapping = V2) %>% mutate(temp = rep(c("22C", "25C"), each = 6, time = 2),nut=rep(c("High", "Low"), each = 12),species = rep(c("BT", "BE", "Tetra","TE","TE","Eup"),time = 4), direc = rep(c("down", "down","single","down","up", "single"), time = 4), dynamics = rep(c("PhenoEco"), time = 24))
PhnEcRho$rho_app <- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_approx_PhnEc ))
PhnEcRho$Lobs_app <- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_approx_PhnEc))

# Pheno-pheno
PhnPhnRho <- str_split_fixed(names(msCCM_approx_PhnPhn), '_', 2) %>% as.data.frame %>% rename(trt = V1,mapping = V2) %>% mutate(temp = rep(c("22C", "25C"), each = 2, time = 2),nut= rep(c("High", "Low"), each = 4),species = rep(c("TE"),time = 8), direc = rep(c("up", "down"), time = 4), dynamics = rep(c("PhenoPheno"), time = 8))
PhnPhnRho$rho_app <- do.call(rbind,mclapply(1:8,Extra.rho,dlist = msCCM_approx_PhnPhn))
PhnPhnRho$Lobs_app <- do.call(rbind,mclapply(1:8,Extra.Lobs,dlist = msCCM_approx_PhnPhn))

rho.all1 <- rbind(EcEcRho,EcPhnRho) %>% full_join(PhnEcRho) %>% full_join(PhnPhnRho) %>% mutate(species = factor(species, levels= c("BT","BE", "TE", "Tetra", "Eup")), direc = factor(direc,levels=c("up", "single","down")),nut = factor(nut,levels = c("Low","High")))

The convergence of cross-mapping skill (i.e. incresing rho with increasing library size) indicates causality. To better understand the causal effects between species interactions, we only used multispatial CCM results that shows convergence to infer those relationships. Here, we removed the CCM results that did not show clear convergence.

# replacing Full25 B->E EcoEco as NA
rho.all1[which(rho.all1$trt=="F25" & rho.all1$mapping=="b.e" & rho.all1$dynamics=="EcoEco"), c("rho_app", "Lobs_app")] <- NA
Appendix II Fig 13
ggplot(data = rho.all1, aes(x=nut, y=rho_app, group=mapping, color= species,shape = direc, lty = direc)) + 
  #geom_hline(yintercept=0,color = "#bcbcbc",  size = 0.3) +
  geom_line(data = rho.all1, aes(x=nut, y=rho_app),size =0.7) +  
  geom_point(size =2, strok =1) +
  ylim(0,1)+
    scale_x_discrete(limits=c("Low","High")) +
    labs(y="CCM Predictability", x = "") +
  facet_wrap(~ dynamics + temp, nrow =2) +
    theme_bw() +
     theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black") ) +
    scale_colour_manual(values=  c("#A997DF","#add0af", "#FFC857","#929084","#2E4052")) +
    scale_linetype_manual(name = "Direction", labels = c("Buttom up", "Single Species","Top down"), values = c(1,3,2)) +
    scale_shape_manual(name = "Direction", labels = c("Buttom up", "Single Species","Top down"), values = c(19, 17, 1))

4.2.2 spline interpolation for MS CCM

4.2.2.1 Multispacial CCM test

As in section 4.2.1, we apply multispacial CCM analysis on time series data with spline interpolation.

# spread the data frame so it's easier to index for CCM
tetra_sp_ccm <- tetra_spline %>%spread( key = treatNum, value = density) %>% arrange(rep,day)
    colnames(tetra_sp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")
eup_sp_ccm<- eup_spline %>%spread( key = treatNum, value = density) %>% arrange(rep,day)
    colnames(eup_sp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")
bac_sp_ccm<- bac_spline %>%spread( key = treatNum, value = density) %>% arrange(rep,day)
    colnames(bac_sp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")

# remove day 15 for tetra and eup to run multispacial CCM with bacteria
    tetra_sp_ccm_d14 <- tetra_sp_ccm[-which(tetra_sp_ccm$day == 15),]
    eup_sp_ccm_d14 <- eup_sp_ccm[-which(eup_sp_ccm$day == 15),]

As before, calculating optimal embeding dimension.

bac_sp_E <- opt_E(6,bac_sp_ccm, "bac")
tetra_sp_E <- opt_E(6,tetra_sp_ccm, "tetra")
eup_sp_E <- opt_E(6,eup_sp_ccm, "eup")
tetra_sp_E_d14 <- opt_E(6,tetra_sp_ccm_d14, "tetra_d14")
eup_sp_E_d14 <- opt_E(6,eup_sp_ccm_d14, "eup_d14")
opt_E_sp1 <- left_join(bac_sp_E,tetra_sp_E, by = "trt") %>% left_join(eup_sp_E, by = "trt") %>% left_join(tetra_sp_E_d14, by = "trt") %>% left_join(eup_sp_E_d14, by = "trt") %>% arrange(trt)

Reformate trait time series that will be used in multispacial CCM analysis.

t.t_sp_ccm <- t.t_sp %>%spread( key = treatNum, value = area) %>% arrange(rep,day) 
    colnames(t.t_sp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")

t.t_sp_ccm_d9 <- t.t_sp_ccm %>% filter(! day %in% c(10:15 ))# only kept tetrahymena trait data from day 0 - day 9 to avoid problems with too many NAs in the time series for ccm test. On average tetrahymena population went extinct by day 9.

e.t_sp_ccm <- e.t_sp %>%spread( key = treatNum, value = area) %>% arrange(rep,day) # only kept tetrahymena trait data from day 0 - day 9 since on average tetrahymena population went extinct by day 9.
    colnames(e.t_sp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")
e.t_sp_ccm_d9 <- e.t_sp_ccm %>% filter(! day %in% c(10:15 ))

# crop ecological dynamics that are interpolated by linear interpolation (approx)
bac_sp_ccm_d9 <- bac_sp_ccm %>% filter(! day %in% c(10:14))
tetra_sp_ccm_d9 <- tetra_sp_ccm %>% filter(! day %in% c(10:14))
eup_sp_ccm_d9 <- eup_sp_ccm %>% filter(! day %in% c(10:14))

Calculate the best embeding dimension for the trait time series, and any cropped time series to match the length of tetrahymena trait time series.

t.t_sp_E <- opt_E(6,t.t_sp_ccm_d9, "t.t")
e.t_sp_E <- opt_E(6,e.t_sp_ccm, "e.t")
e.t.d9_sp_E <- opt_E(6,e.t_sp_ccm_d9, "e.t_d9")
bac.d9_sp_E <- opt_E(6,bac_sp_ccm_d9, "bac_d9")
tetra.d9_sp_E <- opt_E(6,tetra_sp_ccm_d9, "tetra_d9")
eup.d9_sp_E <- opt_E(6,eup_sp_ccm_d9, "eup_d9")
opt_E_sp <- left_join(opt_E_sp1, t.t_sp_E) %>% left_join(e.t_sp_E, by = "trt") %>% left_join(e.t.d9_sp_E, by = "trt") %>% left_join(bac.d9_sp_E, by = "trt") %>% left_join(tetra.d9_sp_E, by = "trt") %>% left_join(eup.d9_sp_E, by = "trt")  %>% arrange(trt)

Here is a list of all the dataframes that will be used for multispacial CCM test with spline interpolation.

dlist_sp <- list(bac_sp_ccm, tetra_sp_ccm,tetra_sp_ccm_d14,eup_sp_ccm,eup_sp_ccm_d14, t.t_sp_ccm, e.t_sp_ccm,e.t_sp_ccm_d9,bac_sp_ccm_d9, tetra_sp_ccm_d9,eup_sp_ccm_d9)
names(dlist_sp) <- c("bac_sp_ccm", "tetra_sp_ccm", "tetra_sp_ccm_d14", "eup_sp_ccm", "eup_sp_ccm_d14", "t.t_sp_ccm", "e.t_sp_ccm","e.t_sp_ccm_d9", "bac_sp_ccm_d9", "tetra_sp_ccm_d9", "eup_sp_ccm_d9")

Eco-Eco interactions. As in section 4.1.1, here we performed the multispacial CCM for the ecological dynamics time series for each of the species pair by treatment.

test_sp_EcEc <- NULL # creating a dataframe to call timeseries for the CCM tests, loop for each treatments
test_sp_EcEc$pair <- c("b.t", "t.b","b.e","e.b","t.e","e.t")
test_sp_EcEc$cause <- c("bac_sp_ccm","tetra_sp_ccm_d14", "bac_sp_ccm","eup_sp_ccm_d14","tetra_sp_ccm","eup_sp_ccm")
test_sp_EcEc$reciever <-  c("tetra_sp_ccm_d14","bac_sp_ccm", "eup_sp_ccm_d14","bac_sp_ccm","eup_sp_ccm","tetra_sp_ccm")
test_sp_EcEc$E <- c("bac","tetra_d14", "bac","eup_d14","tetra","eup")
test_sp_EcEc <- test_sp_EcEc %>% as.data.frame()
msCCM_spline_EcEc <- msCCM_list(dlist_sp, test_sp_EcEc, opt_E_sp, 800)

Eco-Pheno interactions.

test_sp_EcPhn <- NULL # creating a dataframe to call time series for the CCM tests, loop for each treatments
test_sp_EcPhn$pair <- c("b.tt", "b.et","t.tt","t.et","e.tt","e.et")
test_sp_EcPhn$cause <- c("bac_sp_ccm_d9","bac_sp_ccm", "tetra_sp_ccm_d9","tetra_sp_ccm","eup_sp_ccm_d9","eup_sp_ccm")
test_sp_EcPhn$reciever  <- c("t.t_sp_ccm","e.t_sp_ccm", "t.t_sp_ccm","e.t_sp_ccm","t.t_sp_ccm","e.t_sp_ccm")
test_sp_EcPhn$E <- c("bac_d9","bac", "tetra_d9","tetra","eup_d9","eup")
test_sp_EcPhn <- test_sp_EcPhn %>% as.data.frame()               
msCCM_spline_EcPhn <- msCCM_list(dlist_sp, test_sp_EcPhn, opt_E_sp, 800)

Pheno-Eco interactions.

test_sp_PhnEc <- NULL # creating a dataframe to call time series for the CCM tests, loop for each treatments
test_sp_PhnEc$pair <- c("tt.b", "et.b","tt.t","et.t","tt.e","et.e")
test_sp_PhnEc$cause <- c("t.t_sp_ccm","e.t_sp_ccm", "t.t_sp_ccm","e.t_sp_ccm","t.t_sp_ccm","e.t_sp_ccm")
test_sp_PhnEc$reciever  <- c("bac_sp_ccm_d9","bac_sp_ccm", "tetra_sp_ccm_d9","tetra_sp_ccm","eup_sp_ccm_d9","eup_sp_ccm")
test_sp_PhnEc$E <- c("t.t", "e.t", "t.t", "e.t", "t.t", "e.t")
test_sp_PhnEc <- test_sp_PhnEc %>% as.data.frame()                          
msCCM_spline_PhnEc <- msCCM_list(dlist_sp, test_sp_PhnEc, opt_E_sp, 800)

Pheno-Pheno interactions.

test_sp_PhnPhn <- NULL 
test_sp_PhnPhn$pair <- c("tt.et", "et.tt")
test_sp_PhnPhn$cause <- c("t.t_sp_ccm","e.t_sp_ccm_d9")
test_sp_PhnPhn$reciever  <- c("e.t_sp_ccm_d9", "t.t_sp_ccm")
test_sp_PhnPhn$E <- c("t.t", "e.t_d9")
test_sp_PhnPhn <- test_sp_PhnPhn %>% as.data.frame()                          
msCCM_spline_PhnPhn <- msCCM_list(dlist_sp, test_sp_PhnPhn, opt_E_sp, 800)
4.2.2.2 Multispactial CCM plots with Spline interpolation
Appendix II Fig 5
# Eco-Eco msCCM
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
{for(i in c(1:4)){
    plot_msccm(msCCM_spline_EcEc, i*6-5,trt_name[i], EcEc_col,EcEc_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = EcEc_col, lty = EcEc_line, legend =c("B -> T","T -> B", "B -> E", "E -> B", "T -> E","E -> T" ), lwd = 2,  cex = 2)
mtext("Eco-Eco multispacial CCM Prediction (spline)", outer = TRUE, cex = 2)}

Appendix II Fig 6
# Eco-Pheno msCCM 
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
{for(i in c(1:4)){
    plot_msccm(msCCM_spline_EcPhn, i*6-5,trt_name[i],EcPhn_col,EcPhn_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = EcPhn_col, lty = EcPhn_line, legend =c( "B -> Tt", "B -> Et", "T -> Tt","T -> Et", "E -> Tt", "E -> Et"), lwd = 2,  cex = 2)
mtext("Eco-Pheno multispacial CCM Prediction (spline)", outer = TRUE, cex = 2)}

Appendix II Fig 7
# Pheno-Eco msCCM 
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))

par(mfrow = c(2,3),oma = c(0, 2, 4, 0))
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
{for(i in c(1:4)){
    plot_msccm(msCCM_spline_PhnEc, i*6-5,trt_name[i],PhnEc_col, PhnEc_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = PhnEc_col, lty = PhnEc_line, legend =c( "Tt -> B", "Et -> B", "Tt -> T","Et -> T", "Tt -> E", "Et -> E"), lwd = 2,  cex = 2)
mtext("Pheno-Eco multispacial CCM Prediction (spline)", outer = TRUE, cex = 2)}

Appendix II Fig 8
#Pheno-Pheno
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
for(i in c(1:4)){
plot(msCCM_spline_PhnPhn[[2*i-1]][["Lobs"]], msCCM_spline_PhnPhn[[2*i-1]][["rho"]],type="l", col = "#FFC857", lwd = 2, lty = 1, xlim=c(5, 80), ylim=c(0,1), xlab="Library length", ylab="rho",cex.lab=1.5, cex.axis=1.5, cex.names = 1.5, main = trt_name[i] ) # b.t
  lines(msCCM_spline_PhnPhn[[2*i]][["Lobs"]], msCCM_spline_PhnPhn[[2*i]][["rho"]], type="l", col = "#FFC857", lwd = 2, lty = 2) 
}
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = c("#FFC857","#FFC857"), lty = c(1,2), legend =c( "Tt -> Et", "Et -> Tt"), lwd = 2,  cex = 2)
mtext("Pheno-Pheno multispacial CCM Prediction (spline)", outer = TRUE, cex = 2)

4.2.2.3 Plot the rho at largest library size
EcEcRho$rho_sp <- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_spline_EcEc))
EcEcRho$Lobs_sp <- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_spline_EcEc))

EcPhnRho$rho_sp <- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_spline_EcPhn ))
EcPhnRho$Lobs_sp <- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_spline_EcPhn))

PhnEcRho$rho_sp <- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_spline_PhnEc ))
PhnEcRho$Lobs_sp <- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_spline_PhnEc))

PhnPhnRho$rho_sp <- do.call(rbind,mclapply(1:8,Extra.rho,dlist = msCCM_spline_PhnPhn))
PhnPhnRho$Lobs_sp <- do.call(rbind,mclapply(1:8,Extra.Lobs,dlist = msCCM_spline_PhnPhn))

rho.all2 <- rbind(EcEcRho,EcPhnRho) %>% full_join(PhnEcRho) %>% full_join(PhnPhnRho)%>% mutate(species = factor(species, levels= c("BT","BE", "TE", "Tetra", "Eup")), direc = factor(direc,levels=c("up", "single","down")),nut = factor(nut,levels = c("Low","High")))

As mentioned before, the convergence of cross-mapping skill (i.e. increasing rho with increasing library size) indicates causality. To better understand the causal effects between species interactions, we only used multispatial CCM results that shows convergence to infer those relationships. Here, we removed the CCM results that did not show clear convergence. We replace the following CCM results with NAs.

EcoEco Full 25C b.t EcoPheno Half 25C b.et PhenoEco Half 22C et.t
PhenoEco Full 22C et.t PhenoPheno Half 25C et.tt

rho.all2[which(rho.all2$trt=="F25" & rho.all1$mapping=="b.t" & rho.all1$dynamics=="EcoEco"), c("rho_sp", "Lobs_sp")] <- NA
rho.all2[which(rho.all2$trt=="H25" & rho.all1$mapping=="b.et" & rho.all1$dynamics=="EcoPheno"), c("rho_sp", "Lobs_sp")] <- NA
rho.all2[which(rho.all2$trt=="H22" & rho.all1$mapping=="et.t" & rho.all1$dynamics=="PhenoEco"), c("rho_sp", "Lobs_sp")] <- NA
rho.all2[which(rho.all2$trt=="F22" & rho.all1$mapping=="et.t" & rho.all1$dynamics=="PhenoEco"), c("rho_sp", "Lobs_sp")] <- NA
rho.all2[which(rho.all2$trt=="H25" & rho.all1$mapping=="et.tt" & rho.all1$dynamics=="PhenoPheno"), c("rho_sp", "Lobs_sp")] <- NA
Appendix II Fig 14
ggplot(data = rho.all2, aes(x=nut, y=rho_sp, group=mapping, color= species,shape = direc, lty = direc)) + 
  #geom_hline(yintercept=0,color = "#bcbcbc",  size = 0.3) +
  geom_line(data = rho.all2, aes(x=nut, y=rho_sp),size =0.7) +  
  geom_point(size =2, strok =1) +
  ylim(0,1)+
    scale_x_discrete(limits=c("Low","High")) +
    labs(y="CCM Predictability", x = "") +
  facet_wrap(~ dynamics + temp, nrow =2) +
    theme_bw() +
    theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black") ) +
    scale_colour_manual(values=  c("#A997DF","#add0af", "#FFC857","#929084","#2E4052")) +
    scale_linetype_manual(name = "Direction", labels = c("Buttom up", "Single Species","Top down"), values = c(1,3,2)) +
    scale_shape_manual(name = "Direction", labels = c("Buttom up", "Single Species","Top down"), values = c(19, 17, 1))

4.2.3 smooth.spline interpolation for Multispcatial CCM

4.2.3.1 Multispcatial CCM test

Formating ecological time series and calculating optimal embedding dimension for multispacial CCM analysis.

# spread the data frame so it's easier to index for CCM
tetra_smsp_ccm <- tetra_sm_sp %>%spread( key = treatNum, value = density) %>% arrange(rep,day)
    colnames(tetra_smsp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")
eup_smsp_ccm <- eup_sm_sp %>%spread( key = treatNum, value = density) %>% arrange(rep,day)
    colnames(eup_smsp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")
bac_smsp_ccm <- bac_sm_sp %>%spread( key = treatNum, value = density) %>% arrange(rep,day)
    colnames(bac_smsp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")

# remove day 15 for tetra and eup to run multispacial CCM with bacteria
    tetra_smsp_ccm_d14 <- tetra_smsp_ccm[-which(tetra_smsp_ccm$day == 15),]
    eup_smsp_ccm_d14 <- eup_smsp_ccm[-which(eup_smsp_ccm$day == 15),]
bac_smsp_E <- opt_E(6,bac_smsp_ccm, "bac")
tetra_smsp_E <- opt_E(6,tetra_smsp_ccm, "tetra")
eup_smsp_E <- opt_E(6,eup_smsp_ccm, "eup")
tetra_smsp_E_d14 <- opt_E(6,tetra_smsp_ccm_d14, "tetra_d14")
eup_smsp_E_d14 <- opt_E(6,eup_smsp_ccm_d14, "eup_d14")
opt_E_smsp1 <- left_join(bac_smsp_E,tetra_smsp_E, by = "trt") %>% left_join(eup_smsp_E, by = "trt") %>% left_join(tetra_smsp_E_d14, by = "trt") %>% left_join(eup_smsp_E_d14, by = "trt") %>% arrange(trt)
t.t_smsp_ccm <- t.t_sp %>%spread( key = treatNum, value = area) %>% arrange(rep,day) 
    colnames(t.t_smsp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")

t.t_smsp_ccm_d9 <- t.t_smsp_ccm %>% filter(! day %in% c(10:15 ))# only kept tetrahymena trait data from day 0 - day 9 to avoid problems with too many NAs in the time series for ccm test. On average tetrahymena population went extinct by day 9.

e.t_smsp_ccm <- e.t_sp %>%spread( key = treatNum, value = area) %>% arrange(rep,day) # only kept tetrahymena trait data from day 0 - day 9 since on average tetrahymena population went extinct by day 9.
    colnames(e.t_smsp_ccm) <- c("day", "rep","species", "F22", "F25", "H22", "H25")
e.t_smsp_ccm_d9 <- e.t_smsp_ccm %>% filter(! day %in% c(10:15 ))

# crop ecological dynamics that are interpolated by linear interpolation (approx)
bac_smsp_ccm_d9 <- bac_smsp_ccm %>% filter(! day %in% c(10:14))
tetra_smsp_ccm_d9 <- tetra_smsp_ccm %>% filter(! day %in% c(10:14))
eup_smsp_ccm_d9 <- eup_smsp_ccm %>% filter(! day %in% c(10:14))

Calculate the best embeding dimension for the trait time series, and any cropped time series to match the length of tetrahymena trait time series.

t.t_smsp_E <- opt_E(6,t.t_smsp_ccm_d9, "t.t")
e.t_smsp_E <- opt_E(6,e.t_smsp_ccm, "e.t")
e.t.d9_smsp_E <- opt_E(6,e.t_smsp_ccm_d9, "e.t_d9")
bac.d9_smsp_E <- opt_E(6,bac_smsp_ccm_d9, "bac_d9")
tetra.d9_smsp_E <- opt_E(6,tetra_smsp_ccm_d9, "tetra_d9")
eup.d9_smsp_E <- opt_E(6,eup_smsp_ccm_d9, "eup_d9")
opt_E_smsp <- left_join(opt_E_smsp1, t.t_smsp_E) %>% left_join(e.t_smsp_E, by = "trt") %>% left_join(e.t.d9_smsp_E, by = "trt") %>% left_join(bac.d9_smsp_E, by = "trt") %>% left_join(tetra.d9_smsp_E, by = "trt") %>% left_join(eup.d9_smsp_E, by = "trt")  %>% arrange(trt)

Here is a list of all the dataframes that will be used for multispacial CCM test.

dlist_smsp <- list(bac_smsp_ccm, tetra_smsp_ccm,tetra_smsp_ccm_d14,eup_smsp_ccm,eup_smsp_ccm_d14, t.t_smsp_ccm, e.t_smsp_ccm,e.t_smsp_ccm_d9,bac_smsp_ccm_d9, tetra_smsp_ccm_d9,eup_smsp_ccm_d9)
names(dlist_smsp) <- c("bac_smsp_ccm", "tetra_smsp_ccm", "tetra_smsp_ccm_d14", "eup_smsp_ccm", "eup_smsp_ccm_d14", "t.t_smsp_ccm", "e.t_smsp_ccm","e.t_smsp_ccm_d9", "bac_smsp_ccm_d9", "tetra_smsp_ccm_d9", "eup_smsp_ccm_d9")

Eco-Eco interactions.

test_smsp_EcEc <- NULL # creating a dataframe to call timeseries for the CCM tests, loop for each treatments
test_smsp_EcEc$pair <- c("b.t", "t.b","b.e","e.b","t.e","e.t")
test_smsp_EcEc$cause <- c("bac_smsp_ccm","tetra_smsp_ccm_d14", "bac_smsp_ccm","eup_smsp_ccm_d14","tetra_smsp_ccm","eup_smsp_ccm")
test_smsp_EcEc$reciever <-  c("tetra_smsp_ccm_d14","bac_smsp_ccm", "eup_smsp_ccm_d14","bac_smsp_ccm","eup_smsp_ccm","tetra_smsp_ccm")
test_smsp_EcEc$E <- c("bac","tetra_d14", "bac","eup_d14","tetra","eup")
test_smsp_EcEc <- test_smsp_EcEc %>% as.data.frame()
msCCM_smsp_EcEc <- msCCM_list(dlist_smsp, test_smsp_EcEc, opt_E_smsp, 800)

Eco-Pheno interactions.

test_smsp_EcPhn <- NULL # creating a dataframe to call time series for the CCM tests, loop for each treatments
test_smsp_EcPhn$pair <- c("b.tt", "b.et","t.tt","t.et","e.tt","e.et")
test_smsp_EcPhn$cause <- c("bac_smsp_ccm_d9","bac_smsp_ccm", "tetra_smsp_ccm_d9","tetra_smsp_ccm","eup_smsp_ccm_d9","eup_smsp_ccm")
test_smsp_EcPhn$reciever  <- c("t.t_smsp_ccm","e.t_smsp_ccm", "t.t_smsp_ccm","e.t_smsp_ccm","t.t_smsp_ccm","e.t_smsp_ccm")
test_smsp_EcPhn$E <- c("bac_d9","bac", "tetra_d9","tetra","eup_d9","eup")
test_smsp_EcPhn <- test_smsp_EcPhn %>% as.data.frame()               
msCCM_smsp_EcPhn <- msCCM_list(dlist_smsp, test_smsp_EcPhn, opt_E_smsp, 800)

Pheno-Eco interactions.

test_smsp_PhnEc <- NULL # creating a dataframe to call time series for the CCM tests, loop for each treatments
test_smsp_PhnEc$pair <- c("tt.b", "et.b","tt.t","et.t","tt.e","et.e")
test_smsp_PhnEc$cause <- c("t.t_smsp_ccm","e.t_smsp_ccm", "t.t_smsp_ccm","e.t_smsp_ccm","t.t_smsp_ccm","e.t_smsp_ccm")
test_smsp_PhnEc$reciever  <- c("bac_smsp_ccm_d9","bac_smsp_ccm", "tetra_smsp_ccm_d9","tetra_smsp_ccm","eup_smsp_ccm_d9","eup_smsp_ccm")
test_smsp_PhnEc$E <- c("t.t", "e.t", "t.t", "e.t", "t.t", "e.t")
test_smsp_PhnEc <- test_smsp_PhnEc %>% as.data.frame()                          
msCCM_smsp_PhnEc <- msCCM_list(dlist_smsp, test_smsp_PhnEc, opt_E_smsp, 800)

Pheno-Pheno interactions.

test_smsp_PhnPhn <- NULL 
test_smsp_PhnPhn$pair <- c("tt.et", "et.tt")
test_smsp_PhnPhn$cause <- c("t.t_smsp_ccm","e.t_smsp_ccm_d9")
test_smsp_PhnPhn$reciever  <- c("e.t_smsp_ccm_d9", "t.t_smsp_ccm")
test_smsp_PhnPhn$E <- c("t.t", "e.t_d9")
test_smsp_PhnPhn <- test_smsp_PhnPhn %>% as.data.frame()                          
msCCM_smsp_PhnPhn <- msCCM_list(dlist_smsp, test_smsp_PhnPhn, opt_E_smsp, 800)
4.2.3.2 Multispactial CCM plots with Smooth Spline interpolation
Appendix II Fig 9
# Eco-Eco
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
{for(i in c(1:4)){
    plot_msccm(msCCM_smsp_EcEc, i*6-5,trt_name[i], EcEc_col,EcEc_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = EcEc_col, lty = EcEc_line, legend =c("B -> T","T -> B", "B -> E", "E -> B", "T -> E","E -> T" ), lwd = 2,  cex = 2)
mtext("Eco-Eco multispacial CCM Prediction (smooth spline)", outer = TRUE, cex = 2)}

Appendix II Fig 10
# Eco-Pheno msCCM 
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
{for(i in c(1:4)){
    plot_msccm(msCCM_smsp_EcPhn, i*6-5,trt_name[i],EcPhn_col,EcPhn_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = EcPhn_col, lty = EcPhn_line, legend =c( "B -> Tt", "B -> Et", "T -> Tt","T -> Et", "E -> Tt", "E -> Et"), lwd = 2,  cex = 2)
mtext("Eco-Pheno multispacial CCM Prediction (smooth spline)", outer = TRUE, cex = 2)}

###### Appendix II Fig 11

# Pheno-Eco msCCM 
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))

par(mfrow = c(2,3),oma = c(0, 2, 4, 0))
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
{for(i in c(1:4)){
    plot_msccm(msCCM_smsp_PhnEc, i*6-5,trt_name[i],PhnEc_col, PhnEc_line)
  }
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = PhnEc_col, lty = PhnEc_line, legend =c( "Tt -> B", "Et -> B", "Tt -> T","Et -> T", "Tt -> E", "Et -> E"), lwd = 2,  cex = 2)
mtext("Pheno-Eco multispacial CCM Prediction (smooth spline)", outer = TRUE, cex = 2)}

Appendix II Fig 12
#Pheno-Pheno
par(mfrow = c(2,3),oma = c(0, 2, 4, 0), bg = "white")
layout(mat = matrix(c(1,2,5,3,4,0), byrow = T, nrow = 2))
for(i in c(1,3,5,7)){
plot(msCCM_smsp_PhnPhn[[i]][["Lobs"]], msCCM_smsp_PhnPhn[[i]][["rho"]],type="l", col = "#FFC857", lwd = 2, lty = 1, xlim=c(5, 80), ylim=c(0,1), xlab="Library length", ylab="rho",cex.lab=1.5, cex.axis=1.5, cex.names = 1.5 ) # b.t
  lines(msCCM_smsp_PhnPhn[[i+1]][["Lobs"]], msCCM_smsp_PhnPhn[[i+1]][["rho"]], type="l", col = "#FFC857", lwd = 2, lty = 2) 
}
plot(1, type = "n", axes=FALSE, xlab="", ylab="")
legend( "bottomleft", inset=0, title="species pair",col = c("#FFC857","#FFC857"), lty = c(1,2), legend =c( "Tt -> Et", "Et -> Tt"), lwd = 2,  cex = 2)
mtext("Pheno-Pheno multispacial CCM Prediction (smooth spline)", outer = TRUE, cex = 2)

4.2.3.3 Plot the rho at largest library size
EcEcRho$rho_smsp<- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_smsp_EcEc))
EcEcRho$Lobs_smsp<- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_smsp_EcEc))

EcPhnRho$rho_smsp<- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_smsp_EcPhn ))
EcPhnRho$Lobs_smsp<- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_smsp_EcPhn))

PhnEcRho$rho_smsp<- do.call(rbind,mclapply(1:24,Extra.rho,dlist = msCCM_smsp_PhnEc ))
PhnEcRho$Lobs_smsp<- do.call(rbind,mclapply(1:24,Extra.Lobs,dlist = msCCM_smsp_PhnEc))

PhnPhnRho$rho_smsp<- do.call(rbind,mclapply(1:8,Extra.rho,dlist = msCCM_smsp_PhnPhn))
PhnPhnRho$Lobs_smsp<- do.call(rbind,mclapply(1:8,Extra.Lobs,dlist = msCCM_smsp_PhnPhn))

rho.all3 <- rbind(EcEcRho,EcPhnRho) %>% full_join(PhnEcRho) %>% full_join(PhnPhnRho)%>% mutate(species = factor(species, levels= c("BT","BE", "TE", "Tetra", "Eup")), direc = factor(direc,levels=c("up", "single","down")),nut = factor(nut,levels = c("Low","High")))

Again,we removed the following CCM results that did not show clear convergence and replaced them with NAs.

EcoEco Full 22C b.t EcoPheno Half 25C b.et PhenoEco Half 22C et.t PhenoPheno Half 25C et.tt

rho.all3[which(rho.all2$trt=="F22" & rho.all1$mapping=="b.t" & rho.all1$dynamics=="EcoEco"), c("rho_smsp", "Lobs_smsp")] <- NA
rho.all3[which(rho.all2$trt=="H25" & rho.all1$mapping=="b.et" & rho.all1$dynamics=="EcoPheno"), c("rho_smsp", "Lobs_smsp")] <- NA
rho.all3[which(rho.all2$trt=="H22" & rho.all1$mapping=="et.t" & rho.all1$dynamics=="PhenoEco"), c("rho_smsp", "Lobs_smsp")] <- NA
rho.all3[which(rho.all2$trt=="H25" & rho.all1$mapping=="et.tt" & rho.all1$dynamics=="PhenoPheno"), c("rho_smsp", "Lobs_smsp")] <- NA
Appendix II Fig 15
ggplot(data = rho.all3, aes(x=nut, y=rho_smsp, group=mapping, color= species,shape = direc, lty = direc)) + 
  #geom_hline(yintercept=0,color = "#bcbcbc",  size = 0.3) +
  geom_line(data = rho.all3, aes(x=nut, y=rho_smsp),size =0.7) +  
  geom_point(size =2, strok =1) +
  ylim(0,1)+
    scale_x_discrete(limits=c("Low","High")) +
    labs(y="CCM Predictability", x = "") +
  facet_wrap(~ dynamics + temp, nrow =2) +
    theme_bw() +
    theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black") ) +
    scale_colour_manual(values=  c("#A997DF","#add0af", "#FFC857","#929084","#2E4052")) +
    scale_linetype_manual(name = "Direction", labels = c("Buttom up", "Single Species","Top down"), values = c(1,3,2)) +
    scale_shape_manual(name = "Direction", labels = c("Buttom up", "Single Species","Top down"), values = c(19, 17, 1))

4.3 Fig 4 Averged rho plot

We first replaced the CCM results that did no show clear convergence with NAs, as indicated in section 4.2.1.3, 4.2.2.3 and 4.2.3.3.

rho.all3[which(rho.all3$trt=="F25" & rho.all3$mapping=="b.e" & rho.all3$dynamics=="EcoEco"), c("rho_app", "Lobs_app")] <- NA

rho.all3[which(rho.all3$trt=="F25" & rho.all1$mapping=="b.t" & rho.all1$dynamics=="EcoEco"), c("rho_sp", "Lobs_sp")] <- NA
rho.all3[which(rho.all3$trt=="H25" & rho.all1$mapping=="b.et" & rho.all1$dynamics=="EcoPheno"), c("rho_sp", "Lobs_sp")] <- NA
rho.all3[which(rho.all3$trt=="H22" & rho.all1$mapping=="et.t" & rho.all1$dynamics=="PhenoEco"), c("rho_sp", "Lobs_sp")] <- NA
rho.all3[which(rho.all3$trt=="F22" & rho.all1$mapping=="et.t" & rho.all1$dynamics=="PhenoEco"), c("rho_sp", "Lobs_sp")] <- NA
rho.all3[which(rho.all3$trt=="H25" & rho.all1$mapping=="et.tt" & rho.all1$dynamics=="PhenoPheno"), c("rho_sp", "Lobs_sp")] <- NA
# rho.all3 already includes NAs replaced in section 4.2.3.3
rho.all4 <- rho.all3 %>% mutate(rho_mean = rowMeans(rho.all3[, c(8,10,12)], na.rm=TRUE))

Here is the main text figure 4.

ggplot(data = rho.all4, aes(x=nut, y=rho_mean, group=mapping, color= species,shape = direc, lty = direc)) + 
  ylim(0,1)+
  geom_line(data = rho.all4, aes(x=nut, y=rho_mean),size =0.7) +  
  geom_point(size =2, strok =1) +
    scale_x_discrete(limits=c("Low","High")) +
    labs(y="CCM Predictability", x = "") +
  facet_wrap(~ dynamics + temp, nrow =2) +
    theme_bw() +
    theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black") )  +
    scale_colour_manual(values=  c("#A997DF","#add0af", "#FFC857","#929084","#2E4052")) +
    scale_linetype_manual(name = "Direction", labels = c("Bottom up", "Single Species","Top down"), values = c(1,3,2)) +
    scale_shape_manual(name = "Direction", labels = c("Bottom up", "Single Species","Top down"), values = c(19, 17, 1))